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To understand methodological features of the detrended fluctuation analysis (DFA) using a lrigher- 
order polynomial fitting, we establish the direct connection between DFA and Fourier analysis. 

Based on an exact calculation of the single-frequency response of the DFA, the following facts 
are shown analytically: (1) in the analysis of stochastic processes exhibiting a power-law scaling 
of the power spectral density (PSD), S(f) ~ a higher-order detrending in the DFA has no 
adverse effect in the estimation of the DFA scaling exponent a, which satisfies the scaling relation 
a = (/3 + 1)/2; (2) the upper limit of the scaling exponents detectable by the DFA depends on the 
order of polynomial fit used in the DFA, and is bounded by m + 1, where m is the order of the 
polynomial fit; (3) the relation between the time scale in the DFA and the corresponding frequency 
in the PSD are distorted depending on both the order of the DFA and the frequency dependence of 
the PSD. We can improve the scale distortion by introducing the corrected time scale in the DFA 
corresponding to the inverse of the frequency scale in the PSD. In addition, our analytical approach 
makes it possible to characterize variants of the DFA using different types of detrending. As an 
application, properties of the detrending moving average algorithm (DMA) are discussed. 

PACS numbers: 05.40.-a, 02.30.Nw, 02.50.Ey, 05.45.Tp 


I. INTRODUCTION 

Long-range correlations and fractal scaling behavior 
have been observed in a remarkably wide variety of sys¬ 
tems, such as physical 0,1, geophysical M, biologi¬ 
cal Hill, and economic systems jl3-ll4j. In the study 
of the time series data observed in such systems, the 
power spectral analysis is a well-established methodolog¬ 
ical framework BHH- By estimating the slope of the 
log-log scaled power spectral density (PSD), a wide range 
of scaling behavior can be characterized. However, it has 
been pointed out that the power spectral analysis may 
provide spurious detection of scaling behavior caused by 
nonstationarity in time series, such as embedded trends 
and heterogeneous statistical properties fl7l - (2l| . To ob¬ 
tain a more accurate estimate of the scaling exponent, the 
detrended fluctuation analysis (DFA) has been proposed 
Ulllaj and has become a widely used method [lOl[ill .I22i— 
|25| . In this method, local trends in the time series are 
eliminated by least-squares polynomial fitting. If mth- 
order polynomials are employed in the DFA, it is referred 
to as mth-order DFA or DFA?n. In addition, variants 
of the DFA using different types of detrending methods 
have been proposed p6l - [33j . The statistical performance 
and the superiority of the DFAs have been shown by a 
number of numerical studies [31. HU, HH. The effects of 
nonstationarity, nonlinear filters and extreme data loss 
on the DFA sca ling behavior have also been studied sys¬ 
tematically [aura 

However, the mathematical basis of the DFA has not 
been well established. To date, analytical studies on the 
DFA have been mostly limited to the first order case 
[37l - l43j . For instance, Taqqu et al derived the direct 
link between the Hurst exponent H of fractional Brown¬ 


ian motion and the scaling exponent estimated by DFA1 
[37) ]. Furthermore, several researchers have analytically 
studied the relation between DFA1 and power spectral 
analysis to a limited extent [38U4T1 |. However, there are 
few analytical arguments on higher-order DFAs. 

In this work, we attempt to achieve a deeper under¬ 
standing of the methodological features of higher-order 
DFAs by introducing an analytical approach using the 
single-frequency response of the DFA. Based on an ex¬ 
act calculation of the single-frequency response function 
and the assumption of stochastic time series with weak 
trends, the direct connection between higher-order DFAs 
and Fourier analysis can be derived. In previous stud¬ 
ies using numerical experiments, a limitation of the de¬ 
tectable scaling exponent in the DFA has been empir¬ 
ically found |19l 1301 ]. In addition, a deviation of the 
crossover position in the DFA from the corresponding 
frequency in the PSD was also empirically found J34[. 
Our approach can provide clear mathematical reasons 
for these properties and a guiding principle to improve 
the DFA methodology. 

The organization of this paper is as follows. In Sec. El 
we review the DFA method. In Sec. mu we mention 
the basic principle of our approach based on the Fourier 
decomposition. In Sec. IIV1 by considering a single¬ 
frequency component, we derive the single-frequency re¬ 
sponse of DFA. In Sec. El using this derived response, we 
show the direct connection between higher-order DFAs 
and Fourier analysis. In Sec. EH we derive the scaling re¬ 
lation between the scaling exponent estimated by higher- 
order DFAs and that by the PSD. In Sec. IVII1 we derive 
the limitation of the scaling exponent detectable by the 
DFA. In Sec. IVIII1 the distortion between the time scale 
in DFA and frequency scale in PSD are studied. Finally, 
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FIG. 1: (Color online) Decomposition of least squares fitting. 
The time series x(t) (top) is given by the sum of three fre¬ 
quency components (three bottom), x'^(t) +x^ (t) -t-aZ 3 ^ (t). 
Black dashed lines show the least squares line of each time 
series. The least squares line p(t) coincides with p^(t) + 
p (2) (t) +p (3) (f). 


Sec. IIXI provides a summary of our results, and discusses 
possible applications of our method. 


II. DETRENDED FLUCTUATION ANALYSIS 
(DFA) 

Here, we briefly review the DFA [13, EH- The standard 
procedure of the DFA is the following: 1) Starting with 
a time series {a^}^ 1 of length N, the {a^} is integrated 
after subtracting the mean from each data point: 


k -1 

Vk = Yl ( Xi ~ (*)) ’ w 

2—0 

where ( x) denotes the sample mean of {xi). 2) The in¬ 
tegrated time series {yk}^ = i is divided into equal-sized, 
non-overlapping segments of length s. 3) In each seg¬ 
ment, the mean-square-deviation from the least squares 
polynomial fit of order m is calculated. The mean-square- 
deviations are then averaged over all segments and its 
square root F(s), referred to as the fluctuation function, 



FIG. 2: (Color online) Parallel shift of a single frequency 
component y(t). 


is calculated: 


F(s) = 



1 1/2 


l N /s\ frt 


(i) 
■Pk 




(2) 

where [-J is the floor function, andpj^ is the least-squares 
fitting polynomial in the 1th segment. Steps (2) and (3) 
are repeated over multiple time scales (window sizes) to 
characterize the relationship between F(s) and s. A lin¬ 
ear relationship on a log-log plot of F(s) as function of s 
indicates the power-law scaling range, in which the fluc¬ 
tuations can be characterized by a scaling exponent a , 
the slope of the linear relation between log F(s) and log s. 


It is known that the scaling exponent a over large time 
scales is related to the power spectrum exponent (3 by 
a = (0 + 1)/2 (a > 0), where the PSD of {xi} is assumed 
as S(f) ~ / _ ' 9 . Therefore, a white noise time series 
with j3 = 0 is characterized by a = 0.5; a long-range 
correlated time series with 0 < /3 < 1 is indicated by 
0.5 < a < 1; and a long-range anti-correlated time series 
with —1 < /3 < 0 is indicated by 0 < a < 0.5. However, 
the analytical derivation of the scaling relation between 
higher-order DFAs and the PSD has not been reported. 


III. DECOMPOSITION OF LEAST SQUARES 
POLYNOMIAL AND MEAN SQUARE 
DEVIATION 


Before deriving the relation between higher-order DFA 
and Fourier analysis, we mention the basic principle of 
our approach. A similar approach has been used by Will- 
son et al to study the relationship between DFAl and 

psd gam. 

Here, let us consider continuous-time signals, and as¬ 
sume that a time series x{t) in the range [tijlg] is ex¬ 
pressed by the summation of a set of basis functions 
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{* (1) (*)»”• ^ (m) W} [Fig- HI : 

m 

x(t) = ( 3 ) 

i=i 

In this case, if we denote a least squares polynomial fit 
of x(t) and in the range [ti, < 2 ] by p{t) and 

{pb)(i)} ) respectively, the following relation holds (see 


Fig. [H and Appendix A for details): 

m 

P(t) = £p W (f)- (4) 

i=l 

Furthermore, if we consider the square deviation from 
the least-squares polynomial, we obtain 


[ {x(t) - p(x)} 2 dt 

Jt 1 



m „t2 2 m—L m «£ 2 

£/ {s W (i)-p W (ar)} dt + 2^2 £ I { xW W -P (,) w} {s ( j) (*) -p (j) (x)} 

i=l £=1 j=i-\- 


dt. 


(5) 


If the cross terms in Eq. (0 are negligible, the square 
deviation of x(t) can be evaluated by the summation 
of the square deviations of the basis functions {a;M(i)}. 
Note that, in the least-squares fitting, the following al¬ 
ways holds: 

f {x^\t)—p^ l \x)\dt = 0 (i = 1,2, • • ■ , m). (6) 

Jti 

Thus, if values of (t ) — pb) ^ and x ti) (£) _ pU) 

(i ^ j) are distributed independently of each other, the 
expected value of the cross term is zero. 

Based on the above facts, we will consider the fre¬ 
quency components of a time series obtained by Fourier 
transform. In the case of the Fourier decomposition, if 
the phase differences between different frequency compo¬ 
nents are uniformly distributed on [0, 2n], the influence of 
cross-products in Eq. 0 are expected to be negligible. In 
stochastic processes displaying power-law spectra, such 
as fractional Brownian motion and fractional Gaussian 
noise this condition is expected to hold. How¬ 

ever, if a time series includes strong deterministic trends, 
the influence of cross terms in Eq. 0 is not negligible. 
In such a case, our approach is not applicable. 


IV. SINGLE-FREQUENCY RESPONSE 
FUNCTION OF DFA 

To derive the single-frequency response of DFA, we 
consider a single frequency component as a continuous¬ 


time signal, 


x(t) = A cos (2nft + 6 ), (7) 

The integrated signal is given by 


c(t) dr 


A sin (2n ft + 9) 

2^/ 


A sin 9 
27 Tf 


( 8 ) 


Because the constant term in Eq. 0 does not affect the 
fluctuation function F(s) [Eq. 0] in DFA, we neglect 
the second term and study the integrated signal as 


y (t) = sin {2nft + 9 ). 


(9) 


To obtain the fluctuation function F’(s), we have to cal¬ 
culate the mean square deviation from the least squares 
polynomial in each interval [(k — 1 )T, kT], where T is the 
window length, and k denotes the number of the window. 
To simplify the calculation, we shift the interval from 
[(k - 1 )T,kT] to [—T/2,T/2], as shown in Fig. H The 
mean square deviation of the signal y(t ) [Eq. 0] from 
its least squares polynomial can be straightforwardly cal¬ 
culated by the following procedure. The least squares 
polynomial of degree m is obtained by minimizing the 
following function, 






4 


I ({^Cb 1 i O' m\ ) — 


r T / 2 r _A_ 

—T/2 \ 27r/ 


sm 


(2nft + 9 j — ^ a* f > dt, 

i=0 J 


( 10 ) 


where 9 = 9 + TrfT(2k — 1). The coefficients {a,;} of 
the polynomial are determined by solving the following 
equations: 


where i,j = 0,1, ■ • • ,m. After the determination of the 
coefficients {a.^}, the mean square deviation >1> 2 is given 

by 


^({«»}) 

da-j 


( 11 ) 


J 


d> 2 (T,/,A,0) = 


i r T/2 ( a 

Tj_ T/2 \2^f 


sm 


m 

(2nft + 9 j - a,;: 


dt. 


( 12 ) 


i—0 


Moreover, if the phase 9 in Eq. (fl2l) is averaged out, d> 2 
can be approximated by 

1 I* 277 

W{T,f,A) = —J <t> 2 (TJ,A,9)d9, (13) 

_2 

where we refer to $ as the single-frequency response 

function. The analytical formulas of and $ are sum- 

_2 

marized in Appendix B. The plots of $ up to fifth- 
order detrending are shown in Fig. [3] When we ana¬ 
lyze a single-frequency component, the square root of 

_2 

<f> (T = s, /, A) provides the analytical prediction of F(s) 
in the DFA, as shown in Fig. [U However, in the range 
log 10 s < 1.0 in Fig. 3] (a), very small deviations of the es¬ 
timated F(s) by the DFA from the analytical prediction 
by $ are seen. As will be shown in the next subsection, 
these deviations come from the continuous-time approx¬ 
imation of the discrete time series. 




A. Effect of discrete time sampling 

Under the continuous-time approximation, we can an¬ 
alytically derive the single-frequency response function 

_2 

of higher-order DFA. However, in the actual time se¬ 
ries analysis, we have to deal with discrete-time signals. 
Therefore, it is necessary to take the effect of the discrete 
time sampling into account. 

To study this, let us consider discrete time series {a;;}, 

Xi = A cos (2nfi + 9 ), (14) 

where z = 0,1, - - - ,N — 1. The integrated series is given 
by 


k-1 

y k = s T J x i = 

i =0 


A cos {7T f(k — 1) + 9} sm(nfk) 
sin(7r/) 


(15) 
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FIG. 3: (Color online) Single-frequency response functions 
4> 2 (T, /, A) of mth-order DFAs (m = 1,2, •••,5). (a) 

4> 2 (T, / = l, A = 1) versus T. (c) 4> 2 (T = 1,/, A = 1) versus 
/• 


FIG. 4: (Color online) Comparison between the DFA re¬ 
sults (dots) and the analytical predictions by 4>(T = s, f = 
1/100, A = 1) (dashed lines). The F(s) are obtained by an¬ 
alyzing a time series of a single-frequency component, Xi = 
cos (27r*/100) (i = 0,1, • • • , 10' — 1), using mth-order DFAs 
(m= 1,2, •• • ,5). 
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FIG. 5: (Color online) Discretization effect evaluated by 
Eq. (1201) (solid line). Blue dashed lines show the approxi¬ 
mation by exp(7r 2 / 2 /3). 


where k = 1,2, ■ • • , TV. When the window size s is suffi¬ 
ciently long such that the trend in each window is approx¬ 
imated by a constant function, the mean square deviation 
of {yk} around the baseline is given by 

Cc f,A) = ( ie ) 

k= 1 l fc=1 ) 

A 2 

= a ■ 2 , , v (17) 

8sm (7r/) 




login * log ,0 s 


On the other hand, the corresponding quantity in the 
continuous case [Eq. J7])] is given by taking the limit T —> 
oo in Eq. m, 


= lim <S> 2 (T,f,A) 
1 —>-00 

A 2 

8 (?r/) 2 ' 


(18) 

(19) 


Therefore, the discretization effect can be evaluated by 
their ratio 


FIG. 6: (Color online) Comparison of analysis methods, (a) 
An example of numerically generated time series of the frac¬ 
tional Gaussian noise with a = 0.7 after adding piecewise 
linear trends. The piecewise linear trends are described by 
dashed lines. The length of the time series is 10 5 . (b) Fluc¬ 
tuation functions F(s) estimated by mth-order DFA (m = 
1, 2, 3). (c) F(s) estimated by the PSD-based method [Eg. 1271 . 
(c) F(s) estimated by the Fourier-based method [Eg. 1261 . (d) 
Superimposition of (b), (c) and (d). Each point of F(s) rep¬ 
resents the mean value of 100 samples. 


c(/) 


$d (f,A) 

7T 2 / 2 

exp —-— 


(7T/) 2 

sin 2 (7r/) 

for / < 1 
2 


( 20 ) 

( 21 ) 


As shown in Fig. [5l the discretization effect is remark¬ 
ably strong for the high frequency range / > 10“ 1 and 
rapidly decreases as the frequency decreases. Therefore, 
the discretization effect could be negligible, if the high 
frequency components (/ > 10“ x ) are not dominant in 
the time series. 


V. DIRECT CONNECTION BETWEEN DFA 
AND FOURIER ANALYSIS 


with the discrete Fourier transform of real-valued time 
series {x 0l xi, ■ ■ • ,xjv_i}, 


N -1 

X k :=Y / x i e~ 2 ^ kilN , ( 22 ) 

i=0 

where j is the imaginary unit and k = 0,1, • • • , N — 1. 
In this case, the frequency of a single component indexed 
by k is given by fk = k/N. If N is odd, the Fourier 
components with fk are described by 


Based on the above mentioned results, the Fourier 
spectrum of a time series can be converted to the fluc¬ 
tuation function F(s) of DFA. To explain this, we start 


x - fc) = A k cos (2n f k i + 6 k ) for k = 1,2, ••• , N \ 

(23) 




















6 


(b) / 


$f// t 


/// /// 


f//// 

iff// 

■"if/ 

y/fZ 

mil 

-■'if/ 

AmM 

/ + DFA 

/if + DFA 

/ f - PSD 

/ f - PSD 

/ (corrected) 



l°g 10 s log 10 s log 10 s log 10 s 

FIG. 7: (Color online) Comparison between mth-order DFA 
(+) and PSD-based method with (O) or without (A) the 
correction factor c(fk) [Eq. (1201) 1. (a) Analysis of fractional 
Gaussian noise with a = 0.1. (b) Analysis of fractional Gaus¬ 
sian noise with a = 0.3. In the left panels of (a) and (b), 
the correction factor c(fk) is assumed to be c(fk) = 1. Each 
point of F(s ) represents the mean value of 100 samples. 


l°Sio ifJ f) 



FIG. 8: (Color online) Illustration of the single-frequency re- 

_2 

sponse function $ (s, /, A) for a fixed / (solid line) and its 
approximation [Eq. (1301) ] (dashed lines). 


where Ak = 2\X k \/N and 6k = argX^; if N is even, the 
Fourier components are described by 


f x/ = A k cos (27 t f k i + 9 k ) for k = 1,2, • • • , N - 1 

| xi = Ak cos (m) for k = N/2 

(24) 

where Ak = 2\Xk\/N and 6k = arg Xk for k ^ N/2, 

and Ak = \Xk\/N and 0k = 0 for k = N/2. Therefore, 

the time series {.t,} can be expressed by the sum of the 
Fourier components: 

Liv/aj 

x i= ^’ ( 25 ) 

k—0 


where = Xq/ n. 

Using Eqs. m and the fluctuation function of 

DFA can be estimated by 


F(s) 


Liv/ S j L(JV-i)/2J 

tift-t E c (h)* 2 (s,i,f k ,A k ,e k ) 

L / -I i=i fc=i 


(28) 

where [-J is the floor function and 0k = 6k+n f k s[2k — 1). 
In the actual calculation of F(s) through Eq. (l26l) . the 
full information of the phase and amplitude given by the 
Fourier analysis is required. Therefore, we refer to this 
as the Fourier-based method. 

_2 

The single-frequency response function $ [Eq. (fl3l) [ 
also provides an estimate of F(s) as 


F(s) 


L(JV-1)/2J _ 

E c Uk)^ 2 {s,fk,A k ) 

k =1 


1/2 


(27) 


In this case, calculation of F(s) only require the squared 


amplitude given by the PSD. Therefore, we refer to this 
as the PSD-based method. 


As shown in Fig. [6l if we analyze the time series of 
the fractional Gaussian noise embedded in piecewise lin¬ 
ear trends [Fig. [G] (a)], both the Fourier-based method 
and the PSD-based method would be in good agree¬ 
ment with DFA results. However, note that these two 
methods have different characteristics. In the PSD-based 
method [Eq. (|27j) ]. the detrending procedure is not explic¬ 
itly included. That is, the F(s) is only given by locally 
weighted smoothing of the PSD. Therefore, if the trend 
component forms a large contribution to the total power 
of the time series, the PSD-based method cannot repro¬ 
duce the DFA result. On the other hand, in the Fourier- 
based method [Eq. (l26l) [. the detrending procedure is de¬ 
scribed by the relative location between Fourier phases 0 k 
and the local window specified by i (see Fig. [I]). However, 
a strong trend component in the time series could make 
the cross terms in Eq. © large. Therefore, the agree¬ 
ment between the DFA and the Fourier-based method 
is not always excellent. Nevertheless, as we will see, our 
approach serves to provide a deeper understanding of the 
methodological features of the DFA. 


It is also important to note that the discretization ef¬ 
fect evaluated by c(f k ) [Eq. ([20]) ] is remarkably large 
when a strongly anti-correlated time series with a close- 
to-zero value of a is analyzed. As shown in the left panel 
of Fig. [7] (a), if an anti-correlated signal with a = 0.1 is 
analyzed and the c(f k ) is not considered, the difference 
between F(s) of the DFA and that of the PSD-based 
method is clearly seen. The reason for this difference is 
that, in the PSD of such time series, the power contri¬ 
bution of the high frequency components is much larger 
than that of the low frequency components. 


















7 



FIG. 9: (Color online) Illustration of the squared fluctuation 

function F 2 (s) of DFA. The F 2 (s) is given by the weighted 

—2 

sum of single-frequency response functions $ (s, /, A(f)) with 
different frequencies /. Here, the weight is assumed as 
A 2 (f) ~ f~ 0 . (a) P + 1 < 2(m + 1). (b) 0 + 1 > 2(m + 1). 


_ 2 _ 

Base on the asymptotic behavior of $ [Eqs. (051) and 

_2 

(HU)], we approximate $ ( s,f,A ) by 


VI. SCALING RELATION BETWEEN 
HIGHER-ORDER DFA AND PSD 


Using the frequency-component based description of 
the DFA [Eq. (1271) ]. we can study the relation between 

higher-order DFAs and PSD. When s <C 1//Cthe single- 

_2 

frequency response function [Eq. (fldl) ] can be ex¬ 
panded as 




A 2 (7r/) 2m s 2(m+1) 

8(2m + 3)rifcLo(2fc + l) 2 


+ o 


^2(m+2) 


)■ 


where m is the order of the least-squares polynomial fit. 
We have confirmed analytically that Eq. (1281) holds when 
to = 1,2, ■ • • ,7. Thus, we conjecture that Eq. (051) holds 
in general. On the other hand, by taking the limit s —> 


*\s,f,A) 


where 


A 2 (Trf) 2m s 2l ' m+1) 

8(2m + 3)nr=o(2fc + 1 ) 2 

* 

, 8t r 2 / 2 


for s < 

f -> f c 

for s > — 
(30) 


{(2TO+3)ni:L 0 (2fc+i) 2 } 1/(2m+2) 

7r 


(31) 



P 

FIG. 10: (Color online) Estimated scaling exponent a by DFA 
versus the scaling exponent /3 of the analyzed time series. 
Here, numerically generated time series displaying the power- 
law scaling of the PSD, S(f) ~ f~ 0 are analyzed by DFAs. 
The length of the time series is N = 10 5 . The value of a is 
estimated by linear regression of log 10 F(s) versus log 10 s in 
the range 10 2 < s < 10 4 . 


This function is illustrated in Fig. |51 

When the PSD function of a stochastic process is given 
by 


Sx(f) = (32) 


in the range [fi, 1/2], the fluctuation function of the DFA 
can be evaluated by assuming Eq. m and c(/) ~ 1 as 


































F\s) 


, 1/2 


* {s,f,A) 


rfc/s i-tr-i 1 

C$ / f 2m s 2(m+1) (A 2 f~P) df + 

J 0 j fc/s J 


A 2 =S x (f) 


df 


1/2 

fc/s 


df 


= 


/, 


2m—/3+1 


2 ?n — /? + 1 


■cl 2 ) 


/, 


-0-i 


P + 1 

I- 


^ +1 - 


2 0+i 

J+i 


8 


(33) 


(34) 


where we set = 7 r 2m / {8(2m + 3) IlfcLo^ + l) 2 } 

and Cm = ( 87 r 2 )" 1 , and assume // ~ 0 and /3 < 2?n + 1. 
When s> 1 and — 1 < j3 < 2m +1, equation (EH1) results 
in 

F 2 (s)~ s 0+1 4=4> F(s)^s^ 1 . (35) 


Therefore, for — 1 < /3 < 2m + 1, we obtain the scaling 
relation, 


P + 1 


(36) 


where a is the scaling exponent estimated by the DFA 
as F(s) ~ s“. This result shows that the scaling relation 
[Eq. (l36l) ] remains independent of the order of the least- 
squares polynomial fit in the DFA. In other words, the 
higher-order detrending in the DFA has no adverse effect 
in the estimation of the scaling exponent. 


VII. LIMITATION OF DETECTABLE SCALING 
EXPONENT 


In the previous section, the scaling relation between 
DFA and PSD scaling exponents are derived under the 
condition of —1 < p < 2m + 1. This condition implies 
the existence of the limitation of the detectable scaling 
exponent by the DFA. When /3 < —1 and s 0, the 
last term in Eq. m is dominant. This indicates a = 0, 
because of F 2 {s) ~ constant ~ s°. Therefore, the lower 
limit of the detectable scaling exponent a is zero. 

On the other hand, if /3 > 2m + 1, the low frequency 
components near the lower limit of the scaling range, 
fi , are dominant in the contribution to the fluctuation 
function F(s). Therefore, if 2m — /? + 1 takes a negative 
value, and fi is sufficiently small, F 2 (s) given by Eq. (ED 
can be evaluated as 


F 2 (s) 


Cg> [ fc/S f^im+l) (Aq f~P) df 
j Si 



s 2(m+l)^ 


j2m-/3+l 

fc/s \ 

_ 2m — /? + 1 

f=fj 

j2m— /3+1 

fl+e\ 

2 TO — j3 + 1 

f=fij 


g 2(m+l) 

s 2(m+l) 


where 0 < e <C f c /s is chosen such that 

rfa/s rfi+e 

/ / 2 m - / 3 d/« / f 2m ~ P df. 

Jfl+e Jfi 

Thus, the DFA scaling exponent is given by a = m + 1. 
This result shows that the upper limit of the detectable 
scaling exponent a by the DFA is equal to m + 1, which 
depends only on the order of polynomial detrending, 
not on /3. In short, the upper limit of the detectable 
scaling exponent is determined by the power-law expo- 



FIG. 11: (Color online) Estimation of crossover location 
in a first-order autoregressive process [Eq. (144[I ] with a = 
0.9391014. Samples of numerically generated time series 
with length 10 5 were analyzed, (a) Fluctuation functions 
F(s) estimated by mth-order DFA (m = 1,2, ••• ,5). (b) 
F(s) estimated by the PSD-based method [Eq. [27] with 
m = 1,2, ,5. The plots in (a) and (b) were vertically 

shifted for improved visibility, (c) F(s) after adjustment by 
the corrected time scale s = log 10 s — log 10 r, where the log 10 r 
is calculated under the condition of S(f ) ~ / -1 (Table Q]). 
The F(s) estimated by the PSD-based method were vertically 
shifted for improved visibility, (d) The power spectral den¬ 
sity (PSD). The analytical prediction of the crossover point 
is f c = —2.0. Each point represents the mean value of 100 
samples. 


(37) 
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FIG. 12: (Color online) Estimation of crossover location in 
a time series described by the superposition of white noise 
and Brownian noise [Eq. (1491) 1. Samples of numerically gen¬ 
erated time series with length 10 s were analyzed, (a) Fluc¬ 
tuation functions F(s) estimated by mth-order DFA (m = 
1,2, ••• ,5). (b) F(s) estimated by the PSD-based method 
[Eq. [27] with m = 1,2, • • • , 5. The plots in (a) and (b) were 
vertically shifted for improved visibility, (c) F(s) after adjust¬ 
ment by the corrected time scale s = log 10 s — log 10 r, where 
the log 10 r is calculated under the condition of S(f) ~ 
(Table[I]). The F(s) estimated by the PSD-based method were 
vertically shifted for improved visibility, (d) The power spec¬ 
tral density (PSD). The analytical prediction of the crossover 
point is f c = —2.5. Each point represents the mean value of 
100 samples. 


nent of the asymptotic behavior in the left tail of the 
<1>(s,/, A) for s -C 1//. This property can be under¬ 
stood intuitively based on the single-frequency response 

of mth-order DFA. As illustrated in Fig. [9] the left tail 

_2 

of $ (s, /, A) shows asymptotic power-law behavior as 
d> 2 (s, /, A) ~ s 2(m+i) jgq. ([28)1 ] . which corresponds to 
F 2 \s) ~ s 2 ( m +!). As illustrated in Fig. [^a), if f3 + 1 < 

2(m +1), the slope (3 + 1 of log 10 F 2 (s) versus log 10 s can 

_2 

be described by the weighted sum of $ (s, /, A(f)) with 
different frequencies under the condition of A 2 (/) ~ f~d. 
On the other hand, as in Fig. (Hb), if /3 + 1 > 2(m + 1), 
the slope /? + 1 cannot be described by the superposition 
of $>, /, A(f)) under the condition of A 2 (/) ~ f~P. In 
this case, the slope of log 10 F 2 (s) versus log 10 s is deter¬ 
mined by the power-law exponent 2(m+1) in the left tail 

of the $ 2 (s,/, A) for f ~ fi- As shown in Fig. [TU] this 
property is confirmed by numerical experiments. 


TABLE I: Analytical calculations of Eq. Under the con¬ 
ditions shown in this table, the exact results are expressed by 
(a — 7 — In 2)/ In 10 — log 10 s, where 7 is the Euler-Mascheroni 
constant. The values of a are summarized in this table. 



S(f) ~ f 

S(f) - r 1 

s{f) ~ r 2 

DFA1 

91/30 

9/4 

599/420 

DFA2 

737/210 

23/8 

3197/1260 

DFA3 

2407/630 

129/40 

41047/13860 

DFA4 

28097/6930 

139/40 

583507/180180 

DFA5 

382201/90090 

257/70 

621259/180180 


TABLE II: Evaluation of the scale distortion in DFA. The 
values of log 10 r are summarized in this table. 



S(f) ~ f 

s(f) - r 1 

S(f) ~ r 2 

DFA1 

0.268498 

-0.071699 

-0.429475 

DFA2 

0.475305 

0.199735 

0.053075 

DFA3 

0.610419 

0.351738 

0.237321 

DFA4 

0.711943 

0.460312 

0.357587 

DFA5 

0.793605 

0.545620 

0.448582 


VIII. SCALE DISTORTION IN DFA 

It is worth to point out that the relation between the 
time scale in DFA and the corresponding frequency scale 
in PSD is distorted depending on the order of DFA and 
on the PSD shape. In the previous studies based on 
numerical experiments, it has been reported that a devi¬ 
ation of the crossover position observed in the DFA from 
the corresponding frequency in the PSD is enlarged, as 
the order of the DFA increases. However, the relation 
between the scale deviation and the PSD shape has not 
been pointed out. 

A time scale s in DFA and the corresponding fre¬ 
quency / in PSD would be related as s = r/f or 
log 10 s = log 10 r — log 10 /, in which r /1 indicates the 
scale distortion in the DFA. To evaluate r, we calculate 


iogio r = (log 10 /(s)> + log 10 s, 


(38) 


where we define (log 10 f(s)) using the single-frequency 
response function d>” as 


A/2 


—2 


/ (log 10 f)* (s,f,A(f))df 

(login/(*)) = -- —2 -' ( 39 ) 


d> 2 (s,f,A(f)) df 


In Eq. (I5U1) . the amplitude spectrum A(f) is determined 
by the PSD of time series, S(f), as A(f) ~ \JS(f). 
By the relation / ~ 1/s, the asymptotic behavior of 
<f > 2 (s, /, A(f)) as a function of / is given by Eq. (1301) . 
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_ 2 _ 2 

FIG. 13: (Color online) Illustration of a single-frequency response functions <!> (s, /, A). (a) $ for a fixed / in mth-order DFA. 

_ 2 _ 2 

The upper limit of the detectable scaling exponent a is determined by the power-law exponent in the left tail of $ (s). (b) 
for a fixed s. The deviation of the peak position from — log s and the skewed shape result in the scale distortion in the DFA. 


Here, assuming s 1, which means that the peak 

—2 

position of $ as a function of / is located near / = 
1/s ■C 1, we approximate Eq, (1551) as 

f‘ OO 

/ (logio /) (a /, Mf)) df 

(iogio f(s)) « —- p— -• (40) 

/ $ (s,f,A(f))df 

J o 

Under some conditions, equation (HU1) can be calculated 
analytically. For instance, in the case of second-order 
DFA and when S(f) ~ /°, we obtain 

,, ,, u 737 7 + ln27r 

<logi » /(s)) = 2iirhno—mo-- log “> s - (41) 

where 7 « 0.57721566 is the Euler-Mascheroni constant. 
Furthermore, under some other conditions, the analytical 
results of Eq. (®D are expressed by 

(iogio f{s)) = - — ^ nl0 ln27r - log 10 s, (42) 

where a is a constant. The values of a are summarized 
in Table HI 

Based on the calculation of log 10 r [Eq. (1551) ]. we define 
the corrected time scale s in DFA as 

logic 5 = - ( lQ gio f ( s )) = lo gio s - lo Sio ( 43 ) 

In the conditions used in Table [I] the values of log 10 r 
are given by (7 + ln27r — a)/ In 10 (Table HD). Note that 
under the assumption of S(f) ~ and s 1, the 
log 10 r can be approximated by a constant value, but in 
a more general situation, the log 10 r is a function of the 
scale s. As seen in Table m the log 10 r increases, as the 
order of DFA increases, and as the scaling exponent /? 
decreases. In general, as the order of DFA increases, the 

log 10 s are shifted to a smaller scale. This shift can be 

_2 

intuitively understood by the shape of $ ( s,f,A ) for a 

fixed s. That is, as shown in Fig. [3jb) , in the plot of 

_2 

4> (T = 1, /, A = 1) as a function of log frequency, the 
deviation of its peak position from / = 1 /T = 1 and its 


skewed shape are causes for the scale distortion in the 
DFA. 

When we study the crossover phenomena of the scaling 
behavior, it is important to consider the scale distortion 
in DFA. To discuss this, we study two examples display¬ 
ing crossover phenomena. The first example is a time 
series {xi} of a first-order autoregressive process (AR1): 

Xi = aXi- l+Ci, (44) 


where 0 < a < 1 is the constant parameter and is white 
Gaussian noise with zero mean and variance a 2 . In this 
process, the PSD is given by 


P(f) = 


pfc 

APo+PY 


where 


and 



P = 



(45) 


(46) 


(47) 


The frequency dependence of this PSD is nearly constant, 
as it is for white noise, for / < / c , and is asymptotically 
proportional to f ~ 2 for / > f c . 

Here, by setting f c = 10 -2 , er 2 = 1 and a = 
exp(—27r/ c ) = 0.9391014 in the AR 1 process [Eq. (El) ], 
samples of the time series are numerically generated, and 
analyzed by DFA. As shown in Fig. UlT a). the locations 
of crossover points observed in the plots of log 10 F(s) 
versus log 10 s are shifted to the right, as the order of 
DFA increases. Moreover, the location of the crossover 
point in DFA5 is observed at log 10 s ~ 2.6, which shows 
a large deviation from the corresponding crossover point 
log 10 s = log 10 (l// c ) = 2.0. To reduce this deviation, 
we can use the corrected time scales [Eq. (|43|) ]. Here, 
the value of log 10 r calculated by assuming S(f) ~ / -1 , 
shown in Table [ill is used as an approximation of log 10 r 
[Eq. (1551) ]. As shown in Fig. ITIT c). if the corrected time 















11 


scales are used, all curves for different order DFAs are 
almost collapsed into a single curve, and provide a more 
accurate estimate of the crossover point corresponding to 
the analytical prediction. 

The second example is a time series described by the 
superposition of white noise and Brownian noise. A 
Brownian-noise time series (a;) ; } is obtained as the in¬ 
tegral of a white-noise time series. When the variance of 
a white-noise time series is erg, the PSD of the integrated 
series is 




Sb(/) = 


( 27 rfy- 


(48) 


Using this integrated series and another white-noise time 
series {a;( W '*} with variance a 2 , we generate the time series 
{Xi} given by 


(w) (B) 

Xi = x) + x) . 


In this process, the PSD of {x^} is given by 


S(f) =a 2 + 


a 


2 

o 


(27r/r 


(49) 


(50) 


and its crossover frequency f c is 


fc = 


q-p 

2na 


(51) 


The frequency dependence of this PSD is asymptotically 
proportional to f ~ 2 for / < / c , and is nearly constant 
(~ /°) for / > f c . 

Here, by setting f c = 10 —2 5 , cr 2 = 1, and ao = 
2naf c = 0.01986918, samples of the time series are nu¬ 
merically generated and analyzed by DFA. As shown in 
Fig. [H the crossover points are observed to shift as in 
the previous example. In contrast, as shown in Fig. 1121 
if we use the corrected time scale log 10 s, all curves for 
different order DFAs are almost collapsed into a sin¬ 
gle curve, and provide a more accurate estimate of the 
crossover point corresponding to the analytical prediction 
- fogio fc = 2.5. 

The above results demonstrate that the corrected time 
scales log 10 5 based on log 10 f calculated by assuming 
S(f) ~ / -1 can improve the crossover point estima¬ 
tion when the observed scaling exponents lie in the range 
0 < /3 < 2 (or 0.5 < a < 1.5). In more general cases, the 
corrected time scales can be estimated by Eq. m- 


IX. SUMMARY AND DISCUSSION 


FIG. 14: (Color online) Single-frequency response functions 

_2 

<f> (T, /, A) of backward detrending moving average algorithm 
(BDMA) and centered DMA (CDMA) methods, compared 
with mth order DFA (m = 1,2,3). (a) 4 >2 (T, / = 1, A = 1) 
versus T. (c) 4> 2 (T = 1, /, A = 1) versus /. 


PSD, S(f) f &, the higher-order detrending in the 
DFA has no adverse effect on the estimation of the DFA 
scaling exponent a which satisfies the scaling relation 
a = (j3 + l)/2. However, there exists the limitation of 
the detectable scaling exponent a, 0 < a < m+ 1, where 
m is the order of the polynomial fit used in the DFA. This 
fact suggests that, in the case of the first-order DFA, the 
observation of scaling behavior with a ~ 2 may be a spu¬ 
rious detection. Therefore, in such a case, careful testing 
is required. It is also important to note the scale dis¬ 
tortion which is present in DFA and which we suggest 
can be reduced by employing the corrected time scale 
in DFA. As illustrated in Fig. m the above-mentioned 
features are clearly understood by the shape of the single¬ 
frequency response function. 

Our approach based on the single-frequency response 
can provide a more general framework for understanding 
the methodological features of the variants of DFA us¬ 
ing different types of detrending. For instance, we can 
consider the DFA using a simple moving average filter¬ 
ing, known as the detrending moving average algorithm 
(DMA) [13, SH- In the case of the backward DMA 
(BDMA) method, the moving average for a window with 
the size s is defined as 

= ( 52 ) 
j =o 

where {yk} is the integrated time series defined by 
Eq. ([I]). Alternatively, in the case of the centered DMA 
(CDMA) method, the moving average for a window with 
the size s is defined as 


In this work, we have derived the direct connection 
between higher-order DFA and Fourier analysis using 
the single-frequency response of the DFA [Eqs. ESI) and 
<EZD], and shown the usefulness of our approach for un¬ 
derstanding the methodological features of DFA. As an 
important result, it is confirmed that, when we analyze 
stochastic processes displaying a power-law scaling of the 


, LU-i)/2J 

^ s) = - Y, y k +i- ( 53 ) 

i=-L«/2J 

where {yk} is the integrated time series defined by 
Eq. U]). Using these moving averages {y^ }, the fluctua¬ 
tion function F(s) is calculated by the root mean square 
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deviation of the integrated time series around the moving 
average {y{. s ^}. 

To study the methodological properties of the BDMA 
and CDMA methods, we can calculate the single- 


_2 

frequency response function $ (T, /, A) of these meth- 

_2 

ods. These exact expressions of $ (T, /, A) are given by 


^2 frr t _ A 2 {27 r 2 / 2 T 2 - 2n fT sm(2n fT) - cos(2t r/T) + l} 

^ 1G7t 4 /‘ 4 T 2 ' 


r 


(54) 


and by 




A 2 {sin(7r/T) - 7r/T} 2 
87t 4 / 4 T 2 


(55) 


_2 

respectively. Figure [H] shows the $ (T, /, A) of the 
BDMA and CDMA methods, compared with standard 
DFA. The plot of the $ 2 (T, A,/) shows similarity be¬ 
tween CMA and first-order DFA, which indicates that 
both methods have a similar frequency response prop¬ 
erty. As discussed in section IVII1 the power-law tail 
_2 

structure of $ (T, /, A) for s <C 1 // determines the upper 
limit of the detectable scaling exponent by DFA. When 
s <C l//Cequations (lYfl) and (l55l) are expanded as 


9 A I 

4> (T,/,A) = —+0(T 4 ), (56) 


and as 

n A2 2 f2rp4 

$ 2 (r, /, A) = 28 J 8 + 0(T 6 ), (57) 

respectively. Therefore, scaling exponents a estimated 
by these methods are bounded by 1 and 2, respectively. 

Our analytical approach is applicable to higher-order 
DFAs and their variants, and can contribute to a deeper 
understanding of the methodological properties of the 
DFA. Our approach could provide valuable information 
for the improvement of the DFA methodology. 
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Appendix A: Decomposition of least-squares 
polynomial fits 

Consider square-integral functions x(t), x^(t) and 
x^ 2 \t), and assume 

x[t) = x^ ( t ) + x ^ ( t ) (Al) 

in the range [ti, * 2 ] - We denote fcth-order least-squares 
polynomials to x(t), x^(t) and x^(t) by p(t)Cp^(t) 


and p^ 2 \t), respectivly, which are described by 

m 

pit) = ^2 a kt k , (A2) 

(A3) 


k=0 
m 

P (i \t) = {i = 1,2), 

k—0 


where {a*,} and { a£ } {i = 1,2) are constant. In this 
case, the following holds: 


Pit) =P ( 1 ) it) +p [ ' 2 , it). 


(2)7 


(A4) 


The proof of Eq. (IA4I) is the following. The coefficients 
{afc} in Eq. (IA2I) are given by minimization of the sum 
of squared vertical residuals, 

I({ak})= [ {.fit)-Pit )} 2 dt. (A5) 

Jti 


Therefore, {a*,} satisfies 
dl{{a k }) 


da k 

which results in 


= 0 (k = 0,1, • • ■ , m), 


[ {f{t)~p{t)}t k dt = 0. 

Jtx 

In addition, {aj^} and { a satisfy 


(A6) 


(A7) 


rt 2 , 


- p W(t)}t k dt = 0, 

(AS) 

-p {2) {t)}t k dt = 0. 

(A9) 


From Eqs. (IA7I) . (IA8I) and (IA9I) . we can obtaion 

f 2 {p{t)-U l Xt)+pW{t)\}t k dt = 0 

Jti 

TTL 

{ a i - ( a i 1} +a i 2) )} / t l+k dt = 0 (A10) 


i —0 


This equation must always hold true for k = 1,2, ■ • • , to. 
Therefore, we obtain C4 = a,- 1 ^ + a ^ (i = 1, • • • , to), i.e. 
Eq. (IA41) . 


J 
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Appendix B: Analytical formulas of <f> 2 and <f> 2 


1. Zeroth-order detrending (m = 0) 


* 2 ( T,f,A,ff) 
$ 2 (T,f,A) 


A 2 |cos(20) (2sin 2 (7r/T) — nfT sin(27r/T)) + cos(27r/T) + 2n 2 f 2 T 2 — 1} 

16?r 4 / 4 T 2 

A 2 (2t t 2 / 2 T 2 + cos(2t r/T) - l) 

16tt 4 / 4 T 2 


(Bl) 

(B2) 


2. First-order detrending (m = l) 


$ 2 (T,/,A,0) = 32 ^ 6r 4 { ( 3 - 4 ^/ 2 ^ 2 ) cos(27r/T - 20) - 6cos(20) - 4 tt 2 / 2 T 2 cos(20) 

+3 cos(2(t r/T + 0)) - 4tt 2 / 2 T 2 cos(2(t r/T + (9)) + 127r/Tsin(27r/T) 

+67r/Tsin(27r/T - 20) - tt 3 / 3 T 3 sin(27r/T - 20) + 67r/Tsin(2(7r/T + 0)) 

-7r 3 / 3 T 3 sin(2(7r/T + 0)) - 6 - 8tt 2 / 2 T 2 + 4tt 4 / 4 T 4 + (6 - 4tt 2 / 2 T 2 ) cos(2tt/T)} (B3) 
¥ 2 (T,/,A) = 16/ ^ 6y4 (2 tt 4 / 4 T 4 - 47r 2 / 2 T 2 - 3 + (3 - 2^ 2 / 2 T 2 ) caa{2nfT) + 6^/T sin(27r./T)} . (B4) 


3. Second-order detrending (m = 2) 


4> 2 (T,/,A,0) 


$ 2 (tj,a) 


327t 8 / 8 T 6 I 6 ( 7r4 ^ T4 - 24^ 2 / 2 ^ 2 + 15) cos(2 tt/T) + 9Ocos(20) + 47 t 6 / 6 T 6 - 18tt 4 / 4 T 4 

+ (—9t r 4 / 4 T 4 + 78 tt 2 / 2 T 2 - 45) cos(2t r/T - 20) + 24 tt 2 / 2 T 2 cos(20) - 36tt 2 / 2 T 2 - 90 
+6tt 4 / 4 T 4 cos(20) - 45 cos( 2 (t r/T + 0)) + 78 tt 2 / 2 T 2 cos(2(t r/T + 0)) 

-97t 4 / 4 T 4 cos(2(7t/T + 0)) + 1807r/Tsin(27r/T) -48 tt 3 / 3 T 3 sin(27r/T) 

-907r/Tsin(27r/T - 20) + 36 tt 3 / 3 T 3 sin(27r/T - 20) - tt 5 / 5 T 5 sin(27r/T - 20) 
-907r/Tsin(2(7r/T + 0)) + 36 tt 3 / 3 T 3 sin(2(7r/T + 0)) - tt 5 / 5 ? 15 sin(2(7r/T + 0))} (B5) 

16^8/876 ( 3 ^ /4T4 “ 24 ^ 2 / 2T2 + 15 ) cos(2^/T) - 6irfT (4tt 2 / 2 T 2 - 15) sin(27r/T) 
+2 tt 6 / 6 T 6 - 9tt 4 / 4 T 4 - 18tt 2 / 2 T 2 - 45} 


(B6) 
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4. Third-order detrending (m = 3) 


<f> 2 (TJ,A,6) = 32nWfWTS { (-8^ 6 / 6 r 6 + 780 t r 4 / 4 T 4 - 5580^ 2 / 2 r 2 + 3150) cos(2 tt/T) 

- (16tt 6 / 6 T 6 - 540tt 4 / 4 T 4 + 2880tt 2 / 2 T 2 - 1575) cos(2t r/T - 20) - 32 tt 6 / 6 T 6 
—315Ocos(20) - 540tt 2 / 2 T 2 cos(20) - 60 tt 4 / 4 T 4 cos(20) - 120tt 4 / 4 T 4 + 4tt 8 / 8 T 8 
-8?t 6 / 6 T 6 cos(20) + 1575 cos(2(t r/T + 0)) - 2880 tt 2 / 2 T 2 cos(2(t r/T + 0)) 

+540tt 4 / 4 T 4 cos(2(t r/T + 0)) - 16 tt 6 / 6 T 6 cos(2(t r/T + 0)) - 720 tt 2 / 2 T 2 
+63007r/Tsin(27r/T) - 2760 tt 3 / 3 T 3 sin(27r fT) + 120tt 5 / 5 T 5 sin(27r/T) 

+31507r/Tsin(27r/T - 20) - 1560 tt 3 / 3 T 3 sin(27r/T - 20) - tt 7 / 7 T 7 sin(2(7r/T + 0)) 

+120tt 5 / 5 T 5 sin(27r/T - 20) - tt 7 / 7 T 7 sin(27r/T - 20) + 120 tt 5 / 5 T 5 sin(2(7r/T + 0)) 
+31507r/Tsin(2(7r/T + 0)) - 1560 tt 3 / 3 T 3 sin(2(7r/T + 0)) - 3150} (B7) 

¥ 2 (T, /, A) = 167r i^i 0r g {30 tt/T (2tt 4 / 4 T 4 - 467r 2 / 2 T 2 + 105) sm(2nfT) + 2n 8 f 8 T 8 - 16^ 6 / 6 T 6 - 60tt 4 / 4 T 4 

- (4tt 6 / 6 T 6 - 390tt 4 / 4 T 4 + 2790tt 2 / 2 T 2 - 1575) cos(2t r/T) - 360 tt 2 / 2 T 2 - 1575} (B8) 

5. Fourth-order detrending (m = 4) 


$ 2 (T, /, A) = - 19 {27r 10 / 10 T 10 - 257t 8 / 8 T 8 - 150n 6 f 6 T 6 - 1575tt 4 / 4 T 4 - 15750tt 2 / 2 T 2 

167T i ^/ i ^T iU 

+5 (tt 8 / 8 T 8 - 264tt 6 / 6 T 6 + 7245tt 4 / 4 T 4 - 36540tt 2 / 2 T 2 + 19845) cos(2t r/T) 

-30t r/T (4tt 6 / 6 ^ 6 - 287tt 4 / 4 T 4 + 3360tt 2 / 2 T 2 - 6615) sin(27r/T) - 99225} (B9) 


6. Fifth-order detrending (m = 5) 


$ 2 (T, /, A) = - ,, 19 {2tt 12 / 12 T 12 - 367r 10 / 10 T 10 - 315 tt 8 / 8 T 8 - 5040tt 6 / 6 T 6 - 85050tt 4 / 4 T 4 

167T i4 / i4 T i ^ 

-1190700tt 2 / 2 T 2 - 3 (27r 10 / 10 T 10 - 1155tt 8 / 8 T 8 + 81480tt 6 / 6 T 6 - 1417500tt 4 / 4 T 4 

+ 6151950tt 2 / 2 T 2 - 3274425) cos(2t r/T) + 210tt/T (tt 8 / 8 T 8 - 168tt 6 / 6 T 6 

+5724tt 4 / 4 T 4 - 51030tt 2 / 2 T 2 + 93555) sin(27r fT) - 9823275} (BIO) 

7. Sixth-order detrending (m = 6) 


$ 2 (T,/,A) = {2?r 14 / 14 T 14 - 49 t r 12 / 12 T 12 - 588n 10 f w T w - 13230tt 8 / 8 T 8 

-330750?t 6 / 6 T 6 - 7640325tt 4 / 4 T 4 - 137525850tt 2 / 2 T 2 - 42tt/T (87r 10 / 10 T 10 
-2670 tt 8 / 8 T 8 + 203040tt 6 / 6 T 6 - 4916835tt 4 / 4 T 4 + 38045700tt 2 / 2 T 2 - 66891825) sin(27r fT) 
+7 ( 7 r 12 / 12 T 12 - 11047r 10 / 10 T 10 + 162810 tt 8 / 8 T 8 - 6875820tt 6 / 6 T 6 + 95582025tt 4 / 4 T 4 

-381704400tt 2 / 2 T 2 + 200675475) cos(2t r/T) - 1404728325} 


(Bll) 
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8. Seventh-order detrending (m = 7) 


$ 2 (T, /, A) = -_A—-{ - 2t r 16 / 16 T 16 + 64 tt 14 / 14 T 14 + 10087t 12 / 12 T 12 + 302407r lo / 1O T 10 
167T i8 / i8 i 10 l 

+ 1039500tt 8 / 8 T 8 + 34927200tt 6 / 6 T 6 + 1021620600tt 4 / 4 T 4 + 22475653200tt 2 / 2 T 2 
-126?r/T (4 tt 12 / 12 T 12 - 23807r 10 / 10 T 10 + 346170 tt 8 / 8 T 8 - 17901180tt 6 / 6 T 6 
+358107750tt 4 / 4 T 4 - 2541889350tt 2 / 2 T 2 + 4347968625) sm(2nfT) 

+ (8tt 14 / 14 T 14 - 15372?r 12 / 12 T 12 + 41693407r 10 / 10 T 10 - 353617110 tt 8 / 8 T 8 

+11373169500tt 6 / 6 T 6 - 138684996450tt 4 / 4 T 4 + 525368393550tt 2 / 2 T 2 

-273922023375) cos(2tt/T) + 273922023375 j (B12) 
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